Non-diffusive transport in plasma turbulence: a fractional diffusion approach * 
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^ Numerical evidence of non-diffusive transport in three-dimensional, resistive pressure-gradient- 

driven plasma turbulence is presented. It is shown that the probability density function (pdf) of 
test particles' radial displacements is strongly non-Gaussian and exhibits algebraic decaying tails. 
To model these results we propose a macroscopic transport model for the pdf based on the use of 
fractional derivatives in space and time, that incorporate in a unified way space-time non-locality 
(non-Fickian transport), non-Gaussianity, and non-diffusive scaling. The fractional diffusion model 
reproduces the shape, and space-time scaling of the non-Gaussian pdf of turbulent transport calcu- 
lations. The model also reproduces the observed super-diffusive scaling. 



' Recent experimental and theoretical evidence indicates that transport in magnetically confined fusion plasmas devi- 
ates from the standard diffusion paradigm. Typical examples include the confinement time scaling in low confinement 
mode plasmas [1,2], perturbative experiments [3-5], and the non-Gaussianity and long-range correlations of fluctu- 
Q , ations [6]. The standard diffusion paradigm breaks down in these cases because it rests on restrictive assumptions 
■ ^ including locality, Gaussianity, lack of long-range correlations, and linearity. In particular, according to Pick's law, 
^»-^,, the fluxes, which contain the dynamical information of the transport process, are assumed to depend only on local 
' quantities, i.e. the gradients of the fields. Also, at a microscopic level, the diffusion paradigm assumes the existence 
of an underlying un-correlated, Gaussian stochastic process, .i.e. a Brownian random walk. 

The need to develop models that go beyond these restrictive assumptions, is the main motivation of this letter that 
has two connected goals. The first goal is to show numerical evidence of non-diffusive transport in pressure-gradient- 
K*" ■ driven plasma turbulence. We do this by integrating test particles in the E x B field obtained from a nonlinear, 
, three-dimensional turbulence model. Test particle studies of this type have the advantage that incorporate in the 
particle trajectories all the physics of the turbulence model. However, this "microscopic" approach has the limitation 
of being time consuming, and potentially redundant in the sense that it tracks individual, particle orbit information 
that from a statistical point view might be irrelevant. This issue takes us to the second goal of this letter which is 
, to propose and test a macroscopic model describing the statistical properties of transport in pressure-gradient-driven 
plasma turbulence. The proposed model is based on the use of fractional derivative operators which, as it will be 



^ explained below, incorporate in a natural, unified way, non-locality in space and time, non-Gaussianity, and anomalous 



diffusion scaling. 



^ [ The underlying instability in pressure-gradient-driven plasma turbulence is the resistive interchange mode, driven 
^H ^' by the pressure gradient in regions where the magnetic field line curvature is negative. In this system, changes in 
O the pressure gradient trigger instabilities at rational surfaces that locally flatten the pressure profile and increase the 
■ gradient in nearby surfaces. This in turn leads to successive instabilities and intermittent, avalanche-like transport 
', [7], which has been observed to cause anomalous diffusion [8]. This instability is the analog of the Raleigh- Taylor 
' instability, extensively studied in fluids, responsible for the gravity driven overturning of a low density fluid laying 
^ ' below a high density fluid. In magnetically confined plasmas the role of gravity is played by the curvature of the 
■ - - ' magnetic field lines which in a cylindrical geometry is always negative and depends only on the radius. 

The turbulence model that we use, describes the coupled evolution of the electrostatic potential $ and pressure p 
in a cylindrical geometry [7] 
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+ V • V + (Ve)-^ vi $ = -v|$ + ^ + ^vl$ (1) 

(^+V.V+(y,)-^jp=^-^+XxVip + X||V|,, (2) 

where the tilde denotes fluctuating quantities (in time and space), and the angular bracket, (), denotes poloidal and 
toroidal angular (flux surface) average. The magnetic flcld Bq is assumed to be on a cylinder with axis along the z-axis. 
The equilibrium density is no, the ion mass is mj, the averaged radius of curvature of the magnetic field lines is fc, 
and the resistivity is r]. The subindex "±" denotes the direction perpendicular to the cylinder's axis, and the subindex 
"11" denotes the z direction. In both Eqs. (1) and (2) there arc dissipativc terms with characteristic coefficients /i (the 
collisional viscosity) and x± (the coUisional cross-field transport). A parallel dissipation term proportional to X||) is 
also included in the pressure equation. This term can be interpreted as a parallel thermal diffusivity. The evolution 
equation of the flux surface averaged pressure is 

d{p) 19 /~\ „ „1 a / 

' -r{Vrp)=So + D--{r^] (3) 



dr r dr \ I r dr \ dr _ 

It contains a source term, Sq, which is only a function of r. This source of particles and heat is due, for instance, 
to neutral beam heating and fueling. In this case, is essentially determined by the beam deposition profile. In 
the present calculations, we assume a parabolic profile, Sq = 5o [l — (r/a)^]. The model parameters used here are 
fi = 0.2 a^/r/f and xi- = 0.025 a^/r/j, where = o^^q/t] is the resistive time and a the minor radius. The rest of 
the parameters in the model can be reduced to two dimensionless quantities the Lundquist number, which is taken 
to be 5 = 10^, and /3o/2£^ = 0.018, where 0q is the value of /3 at the magnetic axis, and s = a/Ro- The numerical 
calculations were carried out using the KITE code [9] with 363 Fourier components to represent the poloidal and 
toroidal angle dependence for each fluctuating component, and a radial grid resolution of Ar = 7.50 x 10~'*a. 

Having computed the electrostatic potential we study transport by following test particle orbits determined from 
the solutions of the E x B equation of motion 



dr 

rf7 B, 



- = ^V$xBo. (4) 



Since the magnetic field is fixed, the turbulence-induce transport is only due to the fluctuating electrostatic potential. 
This electrostatic approximation, which is quite reasonable for low /? values, is needed in order to carry out the 
numerical calculations in the time range required for reliable transport studies. As an initial condition we used 25, 000 
tracer particles with random initial positions in and z, and radial position r = 0.5 a. Finite size effects did not seem 
to be relevant because during the evolution there were very few particles moving out of the system. In the numerical 
integration of Eq. (4), it is observed that tracer particles either get trapped in eddies for long times, or jump over 
several sets of eddies in a single flight, giving rise to anomalous diffusion [8]. 

Due to the intrinsic stochasticity of test particle orbits, one has to resort to a statistical approach to study transport 
in this system. Our main object of study here is the probability density function (pdf) of radial displacements of the 
particles, P{x,t), where x = (r — a/2)/a and t = t/tji. By definition, at t = 0, P{x,t) = Sx. As t evolves, the pdf 
broadens and develop tails. The triangles in Figure 1 show P{x, t) a.t t = 0.64 obtained from the histogram of particle 
displacements. The log-normal scale of the plot makes clear the strong non-Gaussianity of the density function (in 
this scale a Gaussian is a parabola). The insert in Fig. 1 shows that the tails exhibit algebraic decay with exponent 
equal to 1.75 ± 0.03. The numerical results show that for times above t = 0.1, the moments of the test particles 
displacements exhibits super-diffusive scaling, (a;") ~ t"'', with v = 0.66 ± 0.02. 

In what follows we show that these numerical results can be quantitatively described with a transport model using 
fractional derivative operators in space and time. The generic form of the proposed model is 

cD^P = X [w- aD^ + w+ ,D^] P + A, (5) 

where A is a source, 



r(m-a)^- 1 {x-yr+'- 



^DtP = ^.Tz:^ 9^ r , '''SA-m dy , (6) 

J a 
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are the left and right Riemann-Liouville fractional derivatives respectively, are weighting factors, and m — 1 < 
a < m with m a positive integer. The operator on the left hand side of Eq. (5) is the Caputo fractional derivative in 
time of order < /3 < 1, 

Despite the apparent complexity of their definition, fractional derivatives are natural generalizations of regular deriva- 
tives. In particular, as expected, for a and (3 integers, these operators reduce to regular derivatives, and results of 
regular calculus extend directly to the fractional domain making the analytical study of fractional equations a tractable 
problem. Further information about the definition and basic mathematical properties of these operators can be found 
in Ref. [10]. 

Fractional derivatives are integro-differential operators that incorporate non-locality in space and time in a natural 
way. In particular, the right hand side of Eq. (5) evaluated at a fixed position x takes into account non-local, spatial 
contributions to the flux from all points located to the left (through aD"), and all points located to the right (through 
xDb) of ^ [I-'-]- The constants control the degree of left-right asymmetry in the transport processes. This is 
relevant to fusion plasmas where asymmetric fiuxes are usually observed. The non-locality in time is incorporated in 
the fractional derivative operator on the left hand side of Eq. (5). Here, only the left derivative is used because, by 
causality, transport can only depend on the past history of the system. In addition to the space-time non-locality, 
the fractional diffusion model exhibits non-diffusive scaling of moments. In an infinite domain, the algebraic decaying 
tails of non-Gaussian distributions lead to divergent moments. However, in physical applications (e.g. Ref. [12]) a 
finite- a; cutoff leads to the finite size scaling (cc") ~ f^" , where ly = P/a. Depending on the value of a and /3, transport 
can be super-diffusive (21/ > 1), sub-diffusive {21/ < 1), or diffusive {2u = 1). 

The physics behind the model in Eq. (5) can be further understood from the close connection between transport 
models and the theory of random walks. The standard diffusion model is a macroscopic description of the Brownian 
random walk which assumes that at fixed time intervals t = T, 2T, . . . nT . . . particles at a microscopic level experience 
an un-correlated random displacement, or jump, £n, with probability Vx, where Vx is assumed to have a finite second 
moment. In a similar way, fractional diffusion models can be viewed as macroscopic descriptions of generalized 
Brownian random walk models known as the Continuous Time Random Walk (CTRW) models [12]. In addition to 
the jump probability density Vx, the CTRW model introduces a waiting time probability function Vi- That is, the time 
between jumps, rather than being fixed as in a Brownian walk, it is drawn from a probability function Vt- The different 
types of CTRW processes, and the resulting macroscopic transport models, can be classified based on the characteristic 
waiting time, T, and the characteristic mean-square jump, cr^, being finite or divergent [12]. Based on this, the model 
(5) involving fractional derivatives in space and time can be understood as a general macroscopic description of an 
underlying microscopic stochastic process in which particles exhibit both, jumps without a characteristic spatial scale, 
and waiting times without a characteristic time scale. The space non-locality is a direct consequence of the existence of 
anomalously large jumps (known also as Levy flights) that connect distant regions in space, and the time non-locality 
is due to the history-dependence introduced in the dynamics by the presence of anomalously large waiting times. 

The fractional diffusion model in Eq. (5) is fairly general, and depending on the values of a, /3, and , different 
transport processes can be modeled, including sub-diffusive transport, super-diffusive transport, and asymmetric 
transport. In what follows we show that for the symmetric, super-diffusive transport observed in pressure-gradient- 
driven turbulence: w+ = w~ = —0.5/ cos(7rQ!/2), a = 3/4, f3 = 1/2, and A = 0. To understand this, consider the 
initial value problem of Eq. (5) in an infinite domain, x G (— oo, oo). Setting a = —oo and 6 = cx) in the left and right 
fractional derivative operators, and introducing the Fourier and Laplace transforms 

/oo /"OO 
P{x,t)e''''' dx , P{x,s)= P{x,t)e-''Ut , (9) 
-oo Jo 

Eq. (5) becomes 

(s^ + X |fcr) Pik, s) = sf^-^P{k, 0) , (10) 

where P{k, 0) is the Fourier transform of the initial condition, and we have used the fact that -oo-DJ e**^^ = {ik)" e^''^, 
and xD^ e**^^ = {—ik)°' e'*^^. The test particle transport studies where done by releasing an ensemble of particles at a 
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fixed radius. Based on this, wc consider an initial condition of the form P(x, 0) = S{x). In this case, P{x, t) becomes 
the Green's function or propagator which according to Eq. (10) can be written as 

1 r°° 

P{x, t) = ^J_ e-''^' E0{-x\krt^)dk . (11) 



where 



n=0 



is the Mittag-Lcflcr function [13,10]. As expected, for a = 2 and /3 = 1, P reduces to a Gaussian. For /? = 1 , 
1 < a < 2, P becomes a symmetric Levy stable distribution [14], and for < /? < 1 , a = 2 it reduces to the solution 
of the sub-diffusion fractional equation [12]. Introducing the space-similarity variable, rj = t~^/°'x, the solution can 
be written as 



1 r°° 

Pix, t) = t-^/" K{n) , K{rj) = - / cos{r]z) E0i-xz")dz . 

I" Jo 



(13) 



The solid line in Fig. 1 shows a plot of this solution with fi = 1/2, a = 3/4, and x = 0.09. The agreement with the 
test particles turbulence simulations (triangles) is good. More precisely, using the asymptotic result K(ri) ~ 
for large ry [13], it follows that P{x,to) ~ x~^^~^°'\ for x » t^^", which for a = 3/4 gives a decay exponent equal to 
1.75, a value in very good agreement with the numerical result, 1.75 ± 0.03, shown in the insert in Fig. 1. 

The index /? determines the time-asymptotic scaling properties of P. To show this, we introduce the time-scaling 
variable C = * and write the solution as 

P = |a;|-i C"'^/" K (C"'^/") . (14) 

Using again the large 77, and also the small ry asymptotic behavior of the function K{ri) it follows that P{xo,t) ~ , 
for t <C ]a;o|"^'^, and P(xo, t) ~ t~^, for t ^ |a;o]"/^. This scaling is verified in Fig. 2 that shows the evolution in time 
of P at a fixed position xq- The analytical solution according to Eq. (11), shown with a solid line, exhibits algebraic 
tails in the small t and large t limits, and the expected peak at intermediate times. The circles and the triangles in 
the figure denote the results obtained from the test particles turbulence simulations. The agreement is good, but not 
as sharp as the one in Fig. 2 due to the numerical limitations in the integration of the turbulence model for large 
times. 

As mentioned before, a and /3 determine the scaling properties of the moments of the test particles displacements. 
In particular, (x") ~ t"-" , with v = (3/a. In the present case, a = 3/4 and (3 = 1/2, implies f = 2/3, a value 
in very good agreement with the one obtained from the test particles turbulence simulation, 1/ = 0.66 ± 0.02. The 
super-diffusive scaling implies an anomalous confinement time scaling, tc ~ a"^^. For the case studied here, tc ~ a^^^, 
a reasonable value in the range of the experimentally determined values which typically deviate from the standard- 
diffusion prediction t ^ [2]. 

Summarizing, in this letter we have presented numerical evidence that test particle transport in three-dimensional, 
resistive, pressure-gradient-driven plasma turbulence exhibits non-diffusive transport. In particular, we have shown 
that the pdf of particles is strongly non-Gaussian and exhibits algebraic tails with a decay exponent 1.75 ±0.03. Also, 
the moments of the test particles displacements exhibits supper-diffusive scaling with v = 0.66 ± 0.02. We proposed 
a macroscopic transport model for the pdf based on the use of fractional derivative operators or order a = 3/4 in 
space, and order /? = 1/2 in time. The model incorporates in a natural, unified way, space-time nonlocality (non- 
Fickian transport), non-Gaussianity, and anomalous diffusion scaling. In good agreement with the turbulent transport 
calculations, the pdf in the fractional model decay with exponent 1 + a = 1.75, the pfd scale in time with exponent 
P = 1/2, and the moments scale with exponent ly — /3/a — 2/3 which implies a confinement time scaling t,. ^ a^^^. 
We have focused on symmetric fractional derivatives (i.e. = w~). However, the phenomenology of asymmetric 
operators is quite interesting, and important in fusion plasmas. For example, we have observed that asymmetric 
fractional derivative operators give rise to ballistic- like propagation of pulses. These results indicate that fractional 
diffusion models might be a useful tool to model rapid propagation phenomena in fusion devices. Another area where 
fractional operators looks promising is in the study of the role of non-diffusive transport in the L-H transition. One 
way to approach this problem is to incorporate fractional diffusion operators into reaction-diffusion systems of the 
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type used in L-H transition studies (e.g. Ref. [15]). A first step in this direction was presented in Ref. [16] where it 
was shown that fractional diffusion gives rise to asymmetric, exponential acceleration of fronts. 
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FIG. 1. Non-Gaussian probability density function of test particles in plasma turbulence. The triangles denote the results 
from the histogram of radial displacements according to the test particle, pressure-gradient-driven turbulence model Eqs. (l)-(4). 
The solid line is the analytical solution in Eq. (13) of the symmetric {w'^ = w^) fractional diffusion transport model in Eq. (5) 
with a = 3/4, 13 — 1/2 and x ~ 0.09. The log-log insert shows the algebraic decay of the left (circles) and right (triangles) 
tails. The straight line in the insert is a fit with the predicted decay exponent, 1 + a = 7/4 
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FIG. 2. Time evolution of the probability density function of test particles in pressure-gradient-driven plasma turbulence. 
The circles and the triangles denote the results from the turbulence model Eqs. (l)-(4). The solid line is the analytical solution 
(14) of the symmetric — w^) fractional diffusion transport model in Eq. (5) with a = 3/4, (3 = 1/2 and x = 0.09. In 
agreement with the asymptotic result, the insert shows that the pdf exhibits algebraic tails with exponent equal to f) = 1/2. 



7 



